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Abstract 

The  effects  of  stoichiometry  on  the  atomic  structure  and  the  related  mechanical  properties  of 
boron  carbide  (B4C)  have  been  studied  using  density  functional  theory  and  quantum 
molecular  dynamics  simulations.  Computational  cells  of  boron  carbide  containing  up  to  960 
atoms  and  spanning  compositions  ranging  from  6.7%  to  26.7%  carbon  were  used  to  determine 
the  effects  of  stoichiometry  on  the  atomic  structure,  elastic  properties,  and  stress-strain 
response  as  a  function  of  hydrostatic,  uniaxial,  and  shear  loading  paths.  It  was  found  that 
different  stoichiometries,  as  well  as  variable  atomic  arrangements  within  a  fixed  stoichiometry, 
can  have  a  significant  impact  on  the  yield  stress  of  boron  carbide  when  compressed  uniaxially 
(by  as  much  as  70%  in  some  cases);  the  significantly  reduced  strength  of  boron  carbide  under 
shear  loading  is  also  demonstrated. 

(Some  figures  may  appear  in  colour  only  in  the  online  journal) 


1.  Introduction 

Boron  carbide  (BC),  due  to  its  extreme  hardness,  low  density, 
and  demonstrated  performance,  has  been  used  as  an  armor 
ceramic  for  many  years  [1].  With  nominal  stoichiometry 
B4C,  the  crystal  structure  consists  of  12-atom  icosahedra 
cross-linked  by  3-atom  chains  as  shown  in  figure  1.  Within 
the  structure,  there  is  a  high  degree  of  compositional  variation 
with  configurations  consisting  of  B12  or  BnC  icosahedra 
(among  others)  linked  by  a  variety  of  3 -atom  chains  such 
as  C-C-C  and  C-B-C.  BC  is  generally  regarded  to  have 
R3m  symmetry  [2],  however  this  can  only  be  true  for  a 
subset  of  the  available  atomic  arrangements  since  placement 
of  even  a  single  carbon  atom  within  an  icosahedron  causes 
a  monoclinic  distortion  of  the  rhombohedral  lattice  thereby 
reducing  the  crystalline  symmetry  [3].  Configurations  of 
different  stoichiometry  from  ideal  B4C  (e.g.  B2.75C  or  B5.6C) 
are  referred  to  as  ‘polytypoids’  and  atomic  configurations  with 
ideal  B4C  stoichiometry,  but  different  arrangement  of  atoms 
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within  the  icosahedra  (or  chains),  are  termed  ‘polytypes’. 
Within  each  icosahedron,  there  exist  two  crystallographically 
unique  sites  termed  ‘polar’  and  ‘equatorial’  as  shown  in 
figure  2  and  hereafter,  any  atom  that  specifically  occupies  a 
polar  or  equatorial  site  within  an  icosahedron  will  be  labeled 
with  a  subscript  ‘p’  or  ‘e’  respectively. 

Experimentally,  BC  can  be  produced  by  several  methods 
such  as  reaction  of  boric  oxide  and  carbon  in  an  electric  arc 
furnace  [4]  or  carbothermal  reduction  of  a  boric  acid-citric 
acid  gel  [5].  The  powders  can  consist  of  a  range  of  boron 
to  carbon  ratios  resulting  in  a  complex  phase  diagram  [6], 
an  example  of  which  is  shown  in  figure  3.  Other  phase 
diagrams  have  also  been  constructed  [7-9].  A  new  phase 
equilibrium  diagram  has  just  recently  been  submitted  for 
publication  [10]  which  suggests  that  stoichiometric  B4C  is  a 
line  compound  with  a  monoclinic  structure  that  is  stable  to 
600  K.  In  addition,  it  is  also  suggested  that  a  rhombohedral 
B13C2  solid  solution  phase  is  the  stable  phase  above  600  K 
from  about  10  to  20  atomic  %C.  The  issue  of  a  B13C2  phase 
has  been  debated  for  some  time  after  first  being  identified  by 
Samsonov  et  al  in  1956  [11].  In  a  recent  review  article  [8] 
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Figure  1.  Icosahedral  structure  of  boron  (blue)  carbide  linked  by  a 
3-atom  carbon  (red)  chain.  For  clarity  of  representation,  only  a 
single  chain  is  included  in  the  image. 


Figure  2.  Icosahedron  with  polar  (yellow)  and  equatorial  (blue) 
positions. 

it  has  been  suggested  that  x-ray  diffraction  analysis  of  a 
series  of  boron-rich  materials  indicates  a  distinct  change  in 
the  c  lattice  parameters  at  about  13  atomic  %C,  the  B13C2 
composition.  Further,  in  that  work,  a  single  crystal  of  BC 
was  chemically  analyzed  by  Raman  spectroscopy  and  it  is 
clear  that  even  ‘single  crystals’  of  BC  can  have  significant 
stoichiometric  variation.  McCuistion  et  al  [12]  compiled 
chemical  compositions  of  a  variety  of  BC  powders  also 
indicating  B/C  ratios  from  3.58  to  4.0.  The  exact  chemical 
composition  of  BC  grains  in  two  commercial  bulk  materials 
was  determined  by  Chen  et  al  [13]  using  Electron  Energy  Loss 
Spectroscopy  (EELS)  with  the  following  results:  B/C  ratios  of 


Temperature  *C 


Figure  3.  Boron  carbide  phase  diagram.  Reproduced  with 
permission  from  [6].  Copyright  1983  Elsevier. 


Table  1.  Experimental  bulk  and  shear  moduli  (GPa)  of  boron 
carbide  as  a  function  of  composition. 


Stoichiometry 

%  Carbon 

Bulk  modulus 

Shear  modulus 

B4C 

20.0 

235 

197 

B4.5C 

18.2 

237 

197 

B5.6C 

15.2 

236 

197 

B6.5C 

13.3 

231 

189 

B7.7C 

11.5 

178 

150 

3.81  zb  0.15  and  3.90  ±  0.07.  These  variations  in  composition 
have  an  effect  on  the  observed  mechanical  properties,  as 
shown  in  table  1,  where  the  experimentally  measured  bulk 
and  shear  moduli  for  samples  ranging  from  B4C  (20%  carbon) 
to  B7.7C  (11.5%  carbon)  are  presented  [8].  The  general  trend 
evidenced  by  the  experimental  data  is  a  reduction  in  stiffness 
as  the  boron  concentration  increases.  There  has  been  much 
work  on  determination  of  the  structure  and  properties  of  BC 
samples  of  varying  composition.  Conde  et  al  [14]  determined 
the  hexagonal  lattice  parameters  of  B4C  samples  ranging  from 
10%  to  20%  C  using  glancing  incidence  x-ray  diffraction 
and  similar  experimental  studies  of  structure  as  a  function 
of  stoichiometry  were  conducted  by  Konovalikhin  and 
Ponomarev  [15]  and  Kwei  and  Morosin  [16].  The  application 
of  solid  state  density  functional  theory  (DFT)  [17]  by  Saal 
et  al  [18]  and  Vast  et  al  [19]  provided  theoretically  determined 
structures  and  energies  for  a  range  of  stoichiometries. 

Although  the  structure  and  formation  enthalpy  of  BC  as 
a  function  of  stoichiometry  has  been  well  documented  [18], 
knowledge  of  the  origins  of  the  deformation  and  damage 
mechanisms,  and  their  relation  to  the  macro-mechanical 
properties,  is  critical.  The  discovery  of  shock-induced 
localized  nanoscale  amorphization  [20]  in  hot  pressed  bulk 
BC  and  its  relationship  to  performance  is  unclear.  In  addition, 
the  origin  of  the  dramatic  loss  in  shear  strength  [21,  22]  of 
B4C  in  plate  impact  experiments  has  not  been  elucidated. 
The  elastic  moduli  of  3  BC  polytypes,  all  with  20% 
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carbon  composition,  were  predicted  using  DFT  by  Taylor 
et  al  [23]  and  Aryal  et  al  [24]  under  hydrostatic  and  uniaxial 
load.  Taylor  et  al  showed  that  among  the  Bi2(C-C-C), 
BnCp(C-B-C),  and  BnCe(C-B-C)  polytypes,  there  was  a 
reduction  in  the  bulk  modulus  from  234  to  222  GPa  when 
going  from  a  C-B-C  to  a  C-C-C  chain  at  fixed  20%  C 
composition  and  polytypism  was  shown  to  affect  the  pressure 
evolution  of  the  elastic  moduli  as  well.  Aryal  et  al  [24] 
presented  similar  results  for  the  elastic  moduli  under  uniaxial 
load  and  also  presented  stress-strain  curves  for  two  BC 
polytypes,  both  containing  20%  carbon  and  showed  that  there 
was  a  slight  difference  in  the  stress-strain  response  between 
the  two  polytypes  studied  in  that  work. 

Since  BC  powders  are  most  likely  a  mixture  of  BC 
stoichiometries  [8],  and  given  the  experimentally  observed 
effect  of  such  compositional  variation  on  the  mechanical 
response  of  BC,  it  is  important  that  the  mechanical  properties 
of  a  diversity  of  compositions  be  well  characterized.  This 
information  can  play  a  role  in  identification  of  ‘soft’ 
configurations  that  may  initiate  failure  in  the  BC  structure 
when  impacted  with  high  velocity  projectiles.  In  this  paper,  we 
present  a  comprehensive  survey,  using  DFT,  of  the  mechanical 
properties  of  15  BC  structures  with  stoichiometries  ranging 
from  B2.75C  to  B14C  as  shown  in  table  2.  The  B2.75C 
stoichiometry  (26%  C)  represents  a  carbon-rich  example 
and,  as  indicated  by  the  phase  diagram  in  figure  3,  would 
precipitate  carbon  in  the  form  of  graphite  at  high  processing 
temperatures.  However,  it  was  included  in  the  current  analysis 
purely  as  an  example  of  an  extremely  carbon  rich  limit  in 
order  to  further  elucidate  the  effects  of  stoichiometry  on 
the  mechanical  properties  of  the  icosahedral  structures  under 
consideration.  For  each  stable  configuration,  the  structure, 
elastic  moduli,  and  stress-strain  response  under  several 
loading  paths  using  quantum  molecular  dynamics  (MD) 
simulations  have  been  determined.  We  include  an  analysis 
of  the  stress-strain  response  of  each  structure  under  shear, 
in  addition  to  hydrostatic  and  uniaxial  loading,  since  shear 
has  been  postulated  to  be  a  contributing  factor  to  the 
pressure-induced  amorphization  [25]  phenomenon  that  has 
been  observed  experimentally  in  BC.  The  computational 
approach  adopted  in  this  work  is  described  in  section  2, 
followed  by  presentation  of  the  results  and  discussion  in 
sections  3  and  4  respectively. 

2.  Computational  approach 

2.7.  Crystal  structure  optimization  algorithm 

A  crystal  structure  optimization  program  was  written,  based 
on  the  L-BFGS  [26]  optimization  algorithm,  where  any 
stress  state  (hydrostatic,  uniaxial,  shear)  can  be  imposed. 
Several  software  packages  (such  as  CP2K  [27])  already  offer 
such  implementations,  however,  many  of  those  programs 
rotate  the  input  coordinates  to  a  different  computational 
orientation  making  interpretation  of  the  stress  along  a  specific 
crystallographic  direction,  a  key  component  of  this  work, 
more  difficult  to  monitor  and  interpret.  To  circumvent  this 
difficulty,  a  program  that  does  not  perform  a  rotation  of  the 
input  orientation  was  written  and  used  in  this  study.  At  each 


Table  2.  Boron  carbide  stoichiometries  used  in  this  study.  The 
subscript  e  and  p  labels  denote  equatorial  and  polar  carbons, 
respectively. 


Structure 

Stoichiometry  %  C 

B„Ce(CCC) 

B2.75C 

26.66 

BnCp(CCC) 

B2.75C 

26.66 

Bi2(CCC) 

B4C 

20.00 

B„Ce(CCB) 

b4c 

20.00 

B„CP(CCB) 

b4c 

20.00 

B„Ce(CBC) 

b4c 

20.00 

B„CP(CBC) 

b4c 

20.00 

Bi2(CCB) 

B6.5C 

13.33 

Bi2(CBC) 

B6.5C 

13.33 

B„Ce(BCB) 

B6.5C 

13.33 

B„CP(BCB) 

B6.5C 

13.33 

B„Ce(BBC) 

B6.5C 

13.33 

B„CP(BBC) 

B6.5C 

13.33 

B 12  (BBC) 

Bi4C 

6.66 

Bi2(BCB) 

Bi4C 

6.66 

optimization  step,  the  energy,  forces,  and  stress  tensor  were 
evaluated  using  the  Perdew-Burke-Ernzerhof  [28]  (PBE) 
functional  in  a  double  zeta  valence  plus  polarization  basis 
set  with  a  plane  wave  cutoff  of  800  Ryd  provided  by  the 
CP2K  [27]  program.  It  should  be  noted  that  no  symmetry 
restrictions  were  imposed,  i.e.,  no  constraints  were  applied 
to  enforce  linearity  of  the  3 -atom  chain.  For  each  system,  a 
2x2x2  computational  supercell  containing  120  atoms  was 
used  in  order  to  minimize  size  effects  introduced  when  using 
smaller  computational  cells.  At  each  optimization  step,  the 
stress  tensor  returned  by  CP2K  was  converted  to  cell  vector 
derivatives,  required  by  the  L-BFGS  algorithm  to  update 
the  lattice  vectors,  using  the  transformation  from  stress  to 
cell  vector  gradients  given  by  Doll  [29].  Optimization  was 
considered  converged  when  the  gradient  norm  of  the  cell 
vector  derivatives  was  below  0.0001  atomic  units. 

2.2.  Elastic  constants 

Elastic  constants  are  related  to  the  second  derivative  of  the 
total  energy  with  respect  to  strain,  via 

1  d2E 

Qj  =  ~  —  ~  (l) 

V  dSidsj  0 

where  V  is  the  unit  cell  volume  and  ij  =  1---6  using 
the  compact  Voigt  notation  (1  =  xr,  2  =  yy,  3  =  zz,  4  = 
yz,  5  =  xz,  6  =  xy).  For  this  work,  a  program  was  written 
that  evaluates  the  second  derivatives  in  equation  (1)  via  a 
finite  difference  of  analytic  first  derivatives  of  the  energy  with 
respect  to  strain  (stress  tensor)  provided  by  the  CP2K  code. 
We  have  monitored  the  change  in  the  elastic  constants  as  a 
function  of  hydrostatic  and  uniaxial  load  and  the  required 
stress  corrections  for  elastic  constants,  Cykh  under  non-zero 
load  were  included  using: 

Bijkl  =  Cijkl  T  2  +  $jk&il  T"  SilCfjk  T  SjlCTik 

-  2 SjdCTij)  (2) 

with  Oij  being  an  element  of  the  stress  tensor  and  Bijki 
representing  the  stress  corrected  effective  elastic  constant  (or 
‘Birch  coefficient’)  [30,  31]. 
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Table  3.  Computed  unit  cell  parameters  using  0  K  geometry  optimization  and  MD  simulation  (values  in  parentheses)  at  298  K.  The  shaded 
gray  area  indicates  that  no  stable  configuration  containing  a  linear  3-atom  chain  was  found.  (Lengths  in  Angstroms,  angles  in  degrees, 
volume  in  cubic  Angstroms). 

Structure 

Formula 

%c 

a 

b 

c 

a 

p 

y 

Volume 

Experiment37 

B5.6C 

15.2 

5.19 

5.19 

5.19 

65.18 

65.18 

65.18 

110.02 

B„Ce(CCC) 

B2.75C 

26.66 

5.14(5.15) 

5.21(5.22) 

5.21(5.22) 

64.31(64.28) 

65.19(65.22) 

65.19(65.22) 

109.28(109.93) 

B„CP(CCC) 

B2.75C 

26.66 

5.05(5.06) 

5.21(5.22) 

5.21(5.22) 

64.86(64.82) 

66.05(66.07) 

66.05(66.07) 

108.93(109.66) 

B12(CCC) 

B4C 

20.00 

5.19(5.21) 

5.19(5.21) 

5.19(5.21) 

66.01(66.01) 

66.01(66.01) 

66.01(66.01) 

112.09(112.81) 

B„Ce(CCB) 

b4c 

20.00 

5.16(5.17) 

5.22(5.23) 

5.22(5.23) 

66.18(66.18) 

66.16(66.16) 

66.16(66.16) 

112.72(113.38) 

B„Cp(CCB) 

b4c 

20.00 

5.05(5.06) 

5.23(5.24) 

5.23(5.24) 

66.17(66.14) 

67.33(67.35) 

67.33(67.35) 

112.27(112.90) 

B„Ce(CBC) 

b4c 

20.00 

5.18(5.18) 

5.22(5.22) 

5.22(5.22) 

64.86(64.84) 

64.97(64.97) 

64.97(64.97) 

110.16(110.77) 

B„Cp(CBC) 

b4c 

20.00 

5.07(5.08) 

5.22(5.22) 

5.22(5.22) 

65.24(65.22) 

66.07(66.08) 

66.07(66.08) 

109.75(110.43) 

B 12  (CCB) 

B6.5C 

13.33 

5.16(5.17) 

5.22(5.22) 

5.22(5.22) 

66.79(66.81) 

67.89(67.83) 

67.89(67.83) 

115.29(115.95) 

B 12  (CBC) 

B6.5C 

13.33 

5.20(5.21) 

5.20(5.21) 

5.20(5.21) 

65.83(65.83) 

65.83(65.83) 

65.83(65.83) 

112.10(112.76) 

B„Ce(BCB) 

B6.5C 

13.33 

B„Cp(BCB) 

B6.5C 

13.33 

B„Ce(BBC) 

B6.5C 

13.33 

B„Cp(BBC) 

B6.5C 

13.33 

B 12  (BBC) 

Bi4C 

6.66 

Bi2(BCB) 

Bi4C 

6.66 

2.3.  Quantum  molecular  dynamics  simulations 

Quantum  molecular  dynamics  simulations  of  several  of 
the  structures  were  conducted  under  hydrostatic,  uniaxial, 
and  shear  loading  paths.  For  the  MD  simulations,  the 
computational  cell  size  was  increased  to  a  4  x  4  x  4 
supercell  (960  atoms)  in  order  to  minimize  size  effects.  This  is 
particularly  important  since  large  simulation  cells  are  required 
to  properly  accommodate  large  stresses  and  strains.  For  each 
simulation,  the  atomic  coordinates  were  integrated  using 
the  leap-frog  algorithm  [32]  with  temperature  and  pressure 
controlled  using  algorithms  due  to  Berendsen  [33].  For  each 
MD  trajectory,  atomic  forces  and  stresses  were  computed  by 
CP2K  using  the  PBE  functional  and  basis  set  as  described 
above.  Each  simulation  was  run  for  5000  time  steps  (1  time 
step  =  1  fs),  resulting  in  a  total  simulation  time  of  5  ps  for 
each  system.  In  order  to  determine  the  stress-strain  curves  for 
uniaxial  and  shear  loading,  small  strains  were  applied  in  the 
desired  direction  at  time  step  t  =  0,  and  all  strains  orthogonal 
to  the  initially  applied  strain  were  allowed  to  relax  during  the 
remainder  of  the  simulation.  The  time  averaged  value  of  the 
constrained  stress  tensor  element  was  used  to  generate  the 
stress-strain  curves. 

3.  Results 

3.1.  Structures  at  zero  stress 

The  lattice  parameters  for  each  structure  resulting  from 
optimization  at  0  K  and  MD  simulation  at  298  K  are  presented 
in  table  3.  In  terms  of  the  structure  at  zero  load,  thermal  effects 
are  minimal,  with  slight  expansion  in  the  vector  lengths  and 
minor  variations  in  the  vector  angles,  resulting  from  inclusion 
of  temperature  effects.  The  distortion  of  the  structure  from 
purely  rhombohedral  symmetry  is  clearly  evident  in  many 
of  the  systems  and  within  each  carbon  concentration,  the  a 
lattice  vector  shows  a  larger  contraction  when  the  carbon  atom 


resides  in  the  polar  site.  It  is  noteworthy  that  many  of  the  low 
carbon  content  stoichiometries  produced  structures  that  were 
either  unstable  elastically  (negative  eigenvalue  in  the  elastic 
constant  tensor)  or  converged  to  a  minimum  energy  structure 
with  a  non-linear  3 -atom  chain.  This  seems  to  contradict 
some  of  the  findings  reported  in  earlier  papers  [18],  however 
it  is  not  clear  from  those  publications  if  symmetry  was 
enforced  in  their  calculations  in  order  to  maintain  linearity 
of  the  3 -atom  chain.  The  bending  of  the  3 -atom  chain  for 
these  configurations  was  verified  using  a  plane  wave  basis 
as  implemented  in  the  solid  state  DFT  software  package 
VASP  [34]  in  place  of  the  mixed  Gaussian/plane  wave  basis 
approach  in  CP2K.  Further,  the  bending  of  the  3 -atom  chain 
occurred  in  both  the  2x2x2  and  4x4x4  supercells 
(indicating  that  the  bending  is  not  an  artifact  of  the  size  of 
the  computational  cell  used  in  this  work)  and  was  also  found 
to  occur  when  using  the  local  density  approximation  in  place 
of  the  PBE  functional.  It  is  known  experimentally  that  at  an 
approximately  8%  carbon  content,  boron  begins  to  precipitate 
from  the  lattice  yielding  mixtures  of  boron  carbide  and  pure 
boron  [35].  The  instability  of  the  linear  3-atom  chain  structure 
for  the  many  of  the  boron-rich  compositions  found  in  this 
work  supports  this  finding. 

In  the  double  zeta  basis  set  used  in  this  work,  the 
polar  icosahedral  position  is  energetically  favored  over  the 
equatorial  site  with  the  polar  configurations  being  33.9, 
53.0,  and  35.7  meV/atom  lower  in  energy  than  their 
equatorial  counterparts  for  the  BnC(CCC),  BnC(CCB), 
and  BnC(CBC)  structures  respectively  when  using  the  0  K 
optimized  configurations.  For  the  B13C2  polytypoids,  both 
containing  13.33%  carbon,  the  CBC  arrangement  is  clearly 
favorable  over  the  CCB  arrangement,  which  is  163  meV /atom 
higher  in  energy. 

3.2.  Elastic  constants 

The  zero  stress  elastic  constants  (with  respect  to  rhombo¬ 
hedral  axes),  bulk  modulus,  shear  modulus,  and  Young’s 
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Table  4.  Computed  zero  pressure  elastic  moduli.  All  values  in  GPa.  The  three  values  for  Young’s  modulus  correspond  to  the  values  along 
the  x,  y,  and  z  stress  axes  where  the  z-axis  is  coincident  with  the  rhombohedrally  oriented  [111]  direction. 


Structure 

Formula 

%C 

C11 

C12 

C13 

Cl4 

C33 

C44 

Bulka 

Shear 

Young’s 

Experiment37 

B5.6C 

15.2 

542.8 

130.6 

63.5 

— 

534.5 

164.8 

236.8 

195.6 

460.1 

B„Ce(CCC) 

B2.75C 

26.66 

501.7 

120.0 

65.8 

26.7 

547.7 

190.7 

235.4 

195.4 

467/513/528 

B„Cp(CCC) 

B2.75C 

26.66 

517.9 

133.1 

58.7 

49.1 

544.4 

173.4 

237.4 

189.9 

464/498/531 

Bi2(CCC) 

B4C 

20.00 

486.9 

188.8 

64.9 

14.7 

518.1 

133.6 

221.1 

173.2 

451/451/504 

B„Ce(CCB) 

b4c 

20.00 

457.8 

119.8 

61.5 

36.3 

536.1 

111.6 

221.8 

162.4 

409/438/520 

B„CP(CCB) 

b4c 

20.00 

470.3 

124.6 

53.6 

44.6 

505.7 

132.5 

217.8 

170.9 

417/441/492 

B„Ce(CBC) 

b4c 

20.00 

518.3 

116.8 

65.9 

30.6 

522.5 

159.6 

234.1 

196.5 

481/522/507 

B„CP(CBC) 

b4c 

20.00 

534.2 

120.2 

58.1 

38.0 

525.7 

168.4 

233.6 

199.7 

494/519/514 

Bi2(CCB) 

B6.5C 

13.33 

395.7 

139.4 

82.8 

61.7 

393.7 

96.6 

202.3 

135.6 

291/381/373 

Bi2(CBC) 

B6.5C 

13.33 

531.3 

105.3 

54.2 

-7.95 

528.5 

167.1 

224.4 

201.4 

506/506/519 

a  Computed  using  Voigt-Reuss-Hill  average. 


modulus  for  each  structure  are  presented  in  table  4.  In 
determining  the  elastic  constants,  the  minimum  energy 
structures  found  using  0  K  optimization  (see  table  3)  were 
used  and  temperature  effects  were  not  included.  This  was  done 
to  ensure  that  optimized  structures  with  maximal  symmetry 
were  used  in  computation  of  the  elastic  constant  tensor  in 
order  to  allow  application  of  the  Born  stability  criterion  [36] 
for  the  3 m  point  group  (discussed  below)  to  identify  structural 
instability  under  load.  Inclusion  of  temperature,  in  the  context 
of  MD,  introduces  asymmetry  via  the  random  velocities 
used  to  initiate  the  MD  trajectories.  It  should  be  noted 
that  only  the  Bi2(CCC)  and  Bi2(CBC)  structures  strictly 
adhere  to  R3m  symmetry  and  the  small  distortions  present 
in  the  other  structures  introduce  additional  non-zero,  albeit 
small,  elements  in  the  elastic  constant  tensor.  As  a  result, 
for  ease  of  comparison  across  the  range  of  structures,  only 
the  six  non-zero  Q/’s  that  are  present  for  R3m  symmetry 
are  presented.  The  elastic  constants  are  similar  in  magnitude 
for  all  structures,  however  Bi2(CCB)  shows  a  considerable 
reduction  in  stiffness  compared  to  the  other  systems.  The 
experimental  values  in  table  4  were  measured  by  McClellan 
et  al  using  a  sample  with  stoichiometry  B5.6C  (15.2%  C)  [37]. 
However,  in  their  work,  the  value  of  the  C14  modulus  was 
indeterminate  due  to  the  hexagonal  symmetry  assumed  in 
determination  of  the  elastic  moduli.  The  theoretical  C14 
moduli  are  comparatively  small  and  that  for  the  Bi2(CBC) 
structure  is  negative  compared  to  the  positive  values  obtained 
for  the  rest  of  the  structures.  The  negative  C14  value  for 
Bi2(CBC)  is  consistent  with  the  result  reported  by  Shirai  [38]. 
The  relationship  of  the  values  reported  here  to  the  so-called 
global  minimum  Young’s  modulus  reported  by  McClellan 
et  al  [37]  is  not  clear.  In  addition,  past  research  on  the 
boron-rich  compositions  has  been  inconclusive  on  the  nature 
of  B13C2.  Early  work  by  Samsonov  [11]  suggested  that  it  was 
a  separate  phase  however  this  is  still  an  open  question  in  our 
opinion. 

The  evolution  of  the  elastic  constants  under  hydrostatic 
and  uniaxial  load  for  a  B11C4  polytypoid,  two  B4C  polytypes, 
and  a  B13C2  polytypoid  are  shown  in  figures  4-11.  For 
the  uniaxial  study,  compression  was  applied  along  an  axis 
coincident  with  the  3-atom  chain  ([111]  direction),  which 
is  the  stiffest  elastic  direction  in  the  structure.  All  of  the 
elastic  constants,  which  include  the  stress  corrections  given  in 


Pressure 


Figure  4.  Stress  dependent  elastic  constants  for  BnCp(CCC)  (26% 
C)  under  hydrostatic  load. 

equation  (2),  become  gradually  stiffer  with  pressure;  however, 
the  C44  modulus  decreases  with  load  in  all  cases  except 
that  for  uniaxial  compression  of  the  extremely  carbon  rich 
BnCp(CCC)  (figure  8)  where  it  remains  essentially  constant. 
This  softening  of  C44  has  been  observed  experimentally  [39] 
and  theoretically  [40]  in  alpha  quartz,  which  is  known 
to  undergo  pressure-induced  amorphization  similar  to  that 
observed  in  BC.  This  softening  of  the  C44  shear  modulus  with 
load  may  play  a  role  in  the  sudden  drop  in  shear  strength  of 
shock  loaded  boron  carbide  [21,  22]. 

3.3.  Born  stability  analysis 

We  have  applied  the  Born  stability  criterion  to  identify  stresses 
at  which  stoichiometries  within  the  BC  structure  may  show 
an  elastic  instability.  Born  showed  that  an  expansion  of  the 
internal  energy  of  a  crystal  in  a  power  series  in  the  strain, 
along  with  the  imposition  of  positivity  of  the  energy,  leads  to 
restrictions  on  the  relative  magnitudes  of  the  elastic  constants 
of  a  stable  crystal  [36,  41].  Each  of  the  elastic  constants 
varies  independently  with  stress,  and  at  some  critical  load, 
the  system  may  reach  a  structural  instability.  BC  is  highly 
anisotropic  elastically,  belonging  to  the  crystallographic  space 
group  R3m  with  6  independent  elastic  constants  {Qj}  and 
imposition  of  the  Born  stability  criterion  leads  to  the  following 
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Figure  5.  Stress  dependent  elastic  constants  for  BnCp(CBC)  (20% 
C)  under  hydrostatic  load. 


Pressure 


Figure  6.  Stress  dependent  elastic  constants  for  Bi2(CCC)  (20%  C) 
under  hydrostatic  load. 


Figure  7.  Stress  dependent  elastic  constants  for  Bi2(CBC)  (13%  C) 
under  hydrostatic  load. 


Figure  8.  Stress  dependent  elastic  constants  for  BnCp(CCC)  (26% 
C)  under  uniaxial  load. 


Stress  (GPa) 

Figure  9.  Stress  dependent  elastic  constants  for  BnCp(CBC)  (20% 
C)  under  uniaxial  load. 


Figure  10.  Stress  dependent  elastic  constants  for  Bi2(CCC)  (20% 
C)  under  uniaxial  load. 
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Figure  11.  Stress  dependent  elastic  constants  for  Bi2(CBC) 
(13%  C)  under  uniaxial  load. 


restrictions  on  the  elastic  constants  for  the  BC  structure  (under 
hydrostatic  load): 

*ii-|5i2|>0  (3) 

(#11  +  #12)^33  —  2#13  *  #13  >  0  (4) 

(#11  -  #12)#44  -  2#h  *  #14  >  0  (5) 

where  we  have  used  the  stress  corrected  coefficients  B 
obtained  from  equation  (2)  for  elastic  moduli  at  non-zero 
load,  as  explained  in  [31].  The  general  procedure  is  to 
compute  the  6  elastic  constants  as  a  function  of  load  with 
evaluation  of  equations  (3)-(5)  at  each  point  to  determine  the 
onset  of  the  instability.  Once  the  initial  instability  has  been 
located,  evaluation  of  the  ‘soft  modes’  of  deformation  (atomic 
displacements  corresponding  to  the  instability)  can  also  be 
determined. 

For  most  of  the  structures  (see  table  3)  there  is  a 
monoclinic  distortion  at  zero  load,  consistent  with  conclusions 
reached  by  Huhn  and  Widom  [10].  However,  for  Bi2(CCC) 
and  Bi2(CBC),  there  is  no  reduction  in  symmetry  and  the 
stability  criteria  presented  above  are  strictly  applicable  in 
these  cases.  In  both  cases,  the  first  two  relations  defined 
in  equations  (3)  and  (4)  remained  positive  over  the  applied 
hydrostatic  and  uniaxial  compression  ranges,  however  the 
condition  given  in  equation  (5),  as  shown  in  figure  12,  reaches 
a  zero  value  at  ~67  GPa  for  the  Bi2(CCC)  polytype  after 
hydrostatic  compression  and  ~62  GPa  for  Bi2(CBC)  when 
compressed  uniaxially  along  the  3-atom  chain.  Equation  (5) 
is  violated  before  the  others  due  to  the  decreasing  magnitude 
of  the  C44  elastic  constant  in  each  structure  as  the  load  is 
increased.  Interestingly,  Bi2(CBC),  which  shows  an  elastic 
instability  under  uniaxial  load,  does  not  show  a  critical  point 
over  the  applied  range  under  hydrostatic  load.  However,  the 
hydrostatic  curve  for  this  structure,  as  shown  in  figure  12, 
is  trending  toward  zero  and  extrapolation  of  the  curve 
suggests  an  instability  will  be  reached  at  ~  160  GPa.  Similarly, 
for  the  Bi2(CCC)  uniaxial  curve,  extrapolation  suggests  an 
elastic  instability  will  occur  at  ~  141  GPa.  The  variation  in 
mechanical  response  to  each  particular  loading  pattern  is 


Stress  (GPa) 

Figure  12.  Bom  stability  condition  (see  equation  (5))  for 
Bi2(CCC)  andBi2(CBC). 


exemplary  of  the  changes  that  can  be  induced  by  variation  in 
the  local  bonding  within  the  crystal. 

Although  Bi2(CCC)  and  Bi2(CBC)  show  an  elastic 
instability  as  the  stress  is  increased,  continual  loading 
along  the  same  path,  beyond  the  instability  point,  results 
in  no  discernible  collapse  of  the  structure.  Specifically, 
the  structure  of  the  unit  cell,  for  stresses  beyond  the 
instability,  contains  linear  3-atom  chains  and  symmetric 
icosahedra  although  the  Born  criterion  indicates  that  there  is 
a  lower  energy  configuration  accessible  beyond  the  critical 
stress  under  hydrostatic  or  uniaxial  load.  This  suggests,  at 
least  qualitatively,  that  other  pathways,  possibly  involving 
shear,  are  necessary  in  order  to  access  these  lower  energy, 
lower  symmetry,  configurations.  In  the  current  and  previous 
work  [24,  42]  the  computed  stresses  accommodated  by 
boron  carbide  without  structural  failure  have  been  much 
higher  than  those  suggested  experimentally,  however,  as 
suggested  previously  [24],  the  conditions  of  the  experiments 
may  be  drastically  different  from  the  idealized  models  used 
computationally.  Large  shear  stresses,  which  can  significantly 
lower  phase  transformation  pressures  [43]  may  be  present 
experimentally  and  result  in  much  lower  critical  stresses 
than  those  observed  computationally  along  purely  hydrostatic 
or  uniaxial  paths.  The  effect  of  shear  on  the  mechanical 
properties  of  BC  will  be  discussed  below. 

3.4.  Equations  of  state 

The  hydrostatic  compression  data  for  BnCp(CCC), 
Bi2(CCC),  BnCp(CBC),  and  Bi2(CBC)  resulting  from 
quantum  MD  simulations  at  298  K,  are  shown  in  figure  13. 
The  pressure  response  in  the  absence  of  shear  is  essentially 
equivalent  for  each  structure.  The  resulting  pressure-volume 
data  was  fitted  to  the  third  order  Birch-Murnaghan  equation 
of  state  [44]  and  the  bulk  modulus  and  pressure  derivative  for 
each  structure  are  presented  in  table  5. 

Experimentally,  there  is  a  reduction  in  the  bulk  modulus 
as  the  carbon  concentration  decreases  (see  table  1)  and  this 
trend  is  also  observed  in  the  computed  bulk  moduli. 
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BnCP(CCC)26%C  □ 

BJ2(CCC)-20%  c  o 

BnCp(CBC)-20%  C  A 

B12(CBC)-13%  C  x 


Pressure  (GPa) 

Figure  13.  Pressure-volume  curves  for  structures  within  each 
stoichiometry. 
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Figure  14.  Time  trace  of  the  non-zero  stress  tensor  element  for 
simulations  of  Bi2(CCC)  with  uniaxial  compressions  £3  =  0.08, 
0.180  and  shear  strain  £4  =  0.02. 


Table  5.  Birch-Murnaghan  equation  of  state  for  structures  within 
each  carbon  content. 


Structure 

Formula 

%C 

Bulk  modulus 
(GPa) 

Pressure 

derivative 

BnC(CCC) 

B2.75C 

26.66 

245.5 

1.7 

Bi2(CCC) 

B4C 

20.0 

225.5 

3.3 

B„CP(CBC) 

b4c 

20.0 

226.7 

3.9 

B12(CBC) 

B6.5C 

13.3 

223.9 

2.7 

3.5.  Stress-strain  curves 

Stress-strain  curves  for  each  structure  were  computed  using 
MD  simulations  at  298  K.  A  primary  concern  in  analysis  of 
the  simulation  results  is  the  relatively  short  simulation  time 
of  5  ps  (5000  time  steps)  used  to  integrate  the  trajectories. 
Representative  time  traces  of  the  constrained  stress  tensor 
element,  for  uniaxial  and  shear  strain,  in  Bi2(CCC)  are 
shown  in  figure  14.  The  constrained  tensor  element  (all 
others  are  elements  are  relaxed  to  zero  stress)  has  reached  an 
equilibrium  value  in  all  cases  in  less  than  1000  time  steps 
and  remains  constant  for  the  remainder  of  the  simulation. 
Therefore  the  simulation  time  of  5  ps  is  sufficiently  long  to 
provide  converged  stress-strain  curves  for  the  systems  treated 
in  this  work. 

Stress-strain  curves  for  uniaxial  compression  along  the 
3-atom  chain  axis  for  BnCe(CCC),  BnCp(CCC),  three 
B4C  polytypes,  and  Bi2(CBC)  are  shown  in  figure  15.  The 
maximum  stress  obtained  is  ~140  GPa  in  BnCp(CBC), 
followed  by  Bi2(CBC).  The  carbon-rich  stoichiometries  with 
26%  carbon  have  curves  that  closely  follow  the  others  in  the 
elastic  region,  however  they  reach  a  failure  stress  at  much 
lower  loads,  ~40  GPa,  which  is  about  3.5  times  less  than  that 
of  some  of  the  other  structures.  Although  such  carbon-rich 
compositions  are  not  relevant  according  to  the  phase  diagram, 
this  still  represents  a  dramatic  demonstration  of  the  effect 
of  atomic  structure  and  stoichiometry  on  the  mechanical 
properties.  For  each  material,  the  3 -atom  chain  remains  linear 
up  to  the  maximum  stress  and  then  an  abrupt  bending  of 
the  chain  occurs  resulting  in  a  loss  of  strength.  Snapshots 


Figure  15.  Stress-strain  curves  for  uniaxial  compression  along  axis 
of  3 -atom  chain. 

of  configurations,  extracted  from  the  MD  trajectory  for  the 
BnCp(CCC)  structure  before  and  at  the  critical  stress  are 
shown  in  figure  16.  The  bending  of  the  3-atom  chain  at  the 
failure  point  is  clearly  evident  in  the  structure. 

In  addition  to  uniaxial  compression,  we  have  also 
simulated  the  strength  of  the  structures  under  shear  loading. 
Experimentally,  shear  loading  has  been  identified  as  a  critical 
mechanism  resulting  in  amorphization  of  boron  carbide  [20] 
and  the  shear  strength  of  boron  carbide  has  been  shown  to 
be  significantly  reduced  at  stresses  above  the  Hugoniot  elastic 
limit  [21,  22].  Stress-strain  curves  under  shear,  for  several 
structures,  are  shown  in  figures  17-19.  For  each  system, 
an  £4  shear  strain  was  applied  incrementally,  until  failure, 
and  this  was  done  using  initial  configurations  with  uniaxial 
compression,  <73,  of  0,  10,  20  and  30  GPa  along  the  3-atom 
chain  axis.  In  this  way,  the  shear  strength  of  the  material 
can  be  simulated  as  a  function  of  the  uniaxial  load  on  the 
system.  For  each  material,  the  shear  strength  is  substantially 
less  than  the  uniaxial  compressive  strength  and  as  the  uniaxial 
load  is  increased  there  is  a  reduction  in  the  shear  strength, 
consistent  with  experimental  observation.  The  reduction  in 
shear  strength  is  not  as  marked  as  what  is  seen  in  shock 
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Figure  16.  Snapshots  of  molecular  dynamics  simulation  of  BnCp(CCC)  under  uniaxial  compression  at  20  GPa  (left)  and  at  the  failure 
stress,  40  GPa.  Failure  is  associated  with  bending  of  the  3-atom  chain  (boron  =  blue,  carbon  =  red). 


Figure  17.  Bi2(CCC)  stress-strain  curves  resulting  from  shear 
strain,  s 4 ,  at  several  values  of  uniaxial  stress,  cr3,  along  the  3-atom 
chain  axis.  The  uniaxial  curve  (also  presented  in  figure  15),  in  the 
absence  of  shear,  is  included  for  reference. 


Figure  18.  BnCp(CCB)  stress-strain  curves  resulting  from  shear 
strain,  £4,  at  several  values  of  uniaxial  stress,  (X3,  along  the  3-atom 
chain  axis.  The  uniaxial  curve  (also  presented  in  figure  15),  in  the 
absence  of  shear,  is  included  for  reference. 

experiments.  Shock  loading  is  a  much  more  rapid  process  than 
the  essentially  static  loading  simulations  done  here,  however 
the  reduction  in  shear  strength  trend  is  still  evident  in  the 
simulated  data. 
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Figure  19.  Bi2(CBC)  stress-strain  curves  resulting  from  shear 
strain,  £4,  at  several  values  of  uniaxial  stress,  (T3,  along  the  3-atom 
chain  axis.  The  uniaxial  curve  (also  presented  in  figure  15),  in  the 
absence  of  shear,  is  included  for  reference. 

In  the  work  of  Yan  et  al  [25]  the  effect  of  nonhydro  static 
stress  on  the  elastic  stability  of  BC  at  high  pressures  was 
examined  experimentally.  In  that  work,  there  was  no  evidence 
of  amorphization  as  the  pressure  was  increased.  However, 
as  the  pressure  was  gradually  decreased,  evidence  of  an 
amorphous  phase  was  observed,  occurring  at  a  pressure 
~20  GPa.  As  shown  in  the  shear  strain  curves  presented  in 
figures  17  and  18,  the  yield  strength  under  shear  load  for 
each  structure  is  ~20  GPa.  Although  a  correlation  between 
the  experiment  and  the  computed  value  can  be  inferred,  the 
implication  of  this  finding  is  currently  not  clear  and  will  be 
the  subject  of  future  work. 

4.  Discussion  and  conclusion 

In  conclusion,  atomic  structure  and  stoichiometry  have  a 
marked  effect  on  the  calculated  mechanical  response  of 
boron  carbide.  As  shown  in  table  3,  placement  of  even 
a  single  carbon  atom  within  an  icosahedron  results  in  a 
monoclinic  distortion  of  the  structure  which  reduces  the 
crystal  symmetry  via  contraction  of  the  a  cell  vector  and 
elongation  of  the  other  cell  axes.  The  Bi2(CCB)  structure 
shows  a  considerable  reduction  in  stiffness,  as  evidenced  by 
its  Cn,  C33,  and  C44  elastic  moduli,  that  are  considerably 
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smaller  than  those  of  the  other  structures  and  the  C14  modulus 
of  the  Bi2(CBC)  structure  is  negative  (unlike  the  remaining 
structures)  with  a  value  of  —7.95  GPa.  For  all  structures, 
regardless  of  stoichiometry,  the  C44  elastic  constant  displays 
a  negative  slope  in  response  to  hydrostatic  load  and  as  a 
result  of  this  C44  pressure  softening,  an  elastic  instability 
exists  at  values  of  67  GPa  and  62  GPa  for  the  Bi2(CCC) 
and  Bi2(CBC)  structures  respectively.  The  maximum  yield 
strength  under  uniaxial  compression  is  ~140  GPa,  obtained 
for  the  BnCp(CBC)  structure  (20%  C),  followed  by  the 
Bi2(CBC)  structure  (13%  C)  that  accommodated  uniaxial 
loads  up  to  ~  120  GPa  before  failure.  Within  all  structures,  the 
collapse  of  the  unit  cell  is  associated  with  an  abrupt  bending 
of  the  3-atom  chain  axis  (see  figure  16).  Also  notable  is 
that  all  of  the  materials,  regardless  of  stoichiometry,  exhibit 
pressure  softening  of  the  C44  modulus,  in  contrast  to  all 
the  other  moduli  that  increase  with  pressure.  The  consistent 
softening  of  the  C44  modulus,  regardless  of  stoichiometry, 
indicates  that  this  is  a  feature  of  the  icosahedral/3-atom  chain 
structural  framework  of  the  BC  atomic  structure  which  is 
common  among  all  of  the  materials.  The  softening  of  the  shear 
moduli,  and  the  inclusion  of  shear  strain  which  lowers  the 
yield  stress  considerably,  suggests  possible  atomic  structure 
mechanisms  for  the  experimentally  observed  reduction  in 
shear  strength  during  shock  loading  experiments  and  the 
formation  of  nano- structured  amorphous  regions  observed  in 
ballistic  impact  experiments  at  pressures  of  about  20-25  GPa. 
In  addition,  as  seen  in  figures  17-19,  the  yield  strength  of 
Bi2(CBC)  under  a  shear  strain  (30-40  GPa)  is  about  twice 
that  of  the  other  structures  (10-20  GPa)  suggesting  that  it  may 
be  the  most  stable  under  shear  loading.  The  softening  of  the 
shear  moduli  is  associated  with  the  formation  of  new  bonds 
between  the  unsaturated  central  atom  in  the  3 -atom  chain 
with  equatorial  atoms  in  neighboring  icosahedra.  Uniaxial  and 
hydrostatic  loading  decreases  the  spacing  between  the  central 
chain  and  equatorial  atoms  and  the  formation  of  these  new 
bonds  results  in  an  energetically  more  favorable  configuration. 
Elastic  constants  are  related  to  the  change  in  configuration 
energy  as  a  function  of  displacement  and  the  reduction  in  the 
C44  modulus  is  driven  by  the  formation  of  these  new  bonds 
between  the  chain  and  icosahedra. 

Computational  results  have  suggested  structural  stability 
under  much  higher  loads  than  what  is  observed  experimen¬ 
tally.  It  has  been  suggested  [24]  that  complex  strain  patterns, 
principally  involving  shear,  are  necessary  to  access  the  lower 
energy,  lower  symmetry,  configurations.  Exploration  of  these 
complex  shear  loading  paths  is  possible  using  quantum 
mechanical  potentials,  as  done  in  this  work,  however  the 
exploration  of  the  six  dimensional  strain  space  using  different 
combinations  of  strain  is  computationally  prohibitive.  This 
calls  for  the  development  of  a  classical  potential  applicable 
to  icosahedral  boron  carbide.  The  Reax  [45]  forcefield  may 
serve  as  a  good  functional  form  as  it  can  accommodate  the 
charge  variation  that  occurs  as  a  function  of  geometry  and, 
through  the  use  of  bond  orders,  can  properly  model  carbon 
atoms  which  exist  in  different  hybridization  states  depending 
on  their  presence  in  chains  or  icosahedra.  The  large  amount 
of  data  generated  for  the  structures  and  elastic  properties 


contained  in  this  work  can  serve  as  a  parameterization  set  for 
the  development  of  such  a  classical  model.  Using  a  classical 
model,  much  larger  simulations  of  boron  carbide  can  be 
performed  and  studies  of  the  structural  response  of  boron 
carbide  under  shock  loading  can  be  performed  and  correlated 
with  the  available  experimental  shock  loading  data. 
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